Complex dynamics in coevolution models with ratio-dependent 

functional response 

Per Arne Rikvold 
Center for Materials Research and Technology, 
National High Magnetic Field Laboratory, and Department of Physics, 
Florida State University, Tallahassee, Florida 32306-4350, USA 
Tel.: +1-850-644-6814, Fax.: +1-850-644-6504, E-mail: pnkvold@fsu.edu 

Abstract 

We explore the complex dynamical behavior of two simple predator-prey models of biological 
coevolution that on the ecological level account for interspecific and intraspecific competition, as 
well as adaptive foraging behavior. The underlying individual-based population dynamics are based 
on a ratio-dependent functional response [W.M. Getz, J. theor. Biol. 108, 623 (1984)]. Analytical 
results for fixed-point population sizes in some simple communities are derived and discussed. 
In long kinetic Monte Carlo simulations we find quite robust, approximate 1// noise in species 
diversity and population sizes, as well as power-law distributions for the lifetimes of individual 
species and the durations of periods of relative evolutionary stasis. Adaptive foraging enhances 
coexistence of species and produces a metastable low-diversity phase and a stable high-diversity 
phase. 

Keywords: Complex dynamics; Biological evolution; Coevolution; Predator-prey model; Func- 
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I. INTRODUCTION 



In recent years, there has been a growing reco gnit ion that processes at the ecological 
and evolutionary scales can be strongly linked (l. Hil. 45 . 48]. As a consequence, several 



approaches have been proposed, which model the complex process of coevolution in a fitness 
landscape that changes with the composition of the community, while spanning wide ranges 
of both temporal and taxonomic scales. Early steps in this direction were simulations of 
parapatric and sympatric speciation \^ and the coupled NK model with population dy- 



namics 



(l8l . 1~9| . More recent contributions include the webworld model 



vitn popun 



31, 32 , the 



simple origination/extinction model of Nunes Amaral and Meyer [28|, the speciation model 
of Rossberg et al. (391. t he matching model of Rossberg et al. Hotl, the familv of models in- 
troduced by Yoshida |47j 
versions of the latter [3j 
large individual-based simulations have also been performed of parapatric and sympatric 



the individual-based tangled-nature model jl, 0J]3 and simplified 
as well as network models [3j, |4], |30fl. Recently, 



speciation [13j, [14| and of adaptive radiation [15 



Many of the models mentioned above are deliberately simple, aiming to elucidate uni- 
versal features that are largely independent of the finer details of the ecological interactions 
and the evolutionary mechanisms. While valuable from this point of view, the departure 
of many such models from mechanisms usually included in models of population dynamics 
and evolution has limited their acceptance in the biological communit y A case in point are 
the tangled-nature model 0, E3] and its simplifications 34, 35, H3, 38,|4l|, 49 . Here, we 



therefore introduce a modification of the latter, in which intra- and interspecific competition 
and adaptive foraging are introduced through ratio- dependent functional response functions. 

Our motivation for this study is primarily a desire to understand the extent to which 
the long-time dynamics of complex coevolution models depend on details of the population 
dynamics. By using individual-based models with mutations, we avoid introducing the 
artificial separation between speciations and population dynamics inherent in all species- 
level models, including those mentioned above Q H B H El, m, B EJ |47j. Despite being 
individual-based, our models enable fast simulation of large communities over time scales of 
tens of millions of generations. An important question that one would like to answer in the 
future is whether the avalanche- like mass extinctions observed in the fossil record are due to 
intrinsic fluctuations of the nonlinear dynamics [so-called self-organized criticality or SOC 
26[], or to external perturbations such as asteroid impacts, volcanic eruptions, or climate 
changes, or to a combination of intrinsic and extrinsic causes 27] . In order to successfully 
address this question, it is necessary to understand better the influence that the population 
dynamics have on the intrinsic fluctuations in well-defined model systems. We find that the 
dynamics of the functional-response models studied in this paper differ from the tangled- 
nature type models studied earlier and also from each other, depending on whether or not 
adaptive foraging is implemented. 

The remainder of this paper is organized as follows. The model without adaptive foraging 
is defined in Sec. [TTJ For this model, analytical results for simple communities are derived 
and discussed in Sec. IIII[ and kinetic Monte Carlo simulations of multispecies communities 
with mutations are performed and analyzed in Sec. [IVl Adaptive foraging is introduced 
and investigated in Sec. [V] both by numerical solution of the steady-state equations for 
two-species communities, and by long-time kinetic Monte Carlo simulations for evolving 
multispecies communities. A summary and conclusions are given in Sec. I VI I 



2 



II. MODELS 



We recently performed detailed analytical and simulational studies of the long-time dy- 
namics and community structures of simplified tangled- nature models 34, 35, 



In particular, the behaviors of mutualistic and predator-prey versions were compared by 
Rikvold j35| . and community structures of the latter were compared with data from real 



ecosystems by Rikvold and Sevim 37]. Here we first describe features of these models that 



are shared by the new models that will be introduced below. More detailed descriptions of 



our previous models are given by Rikvold [35j and Rikvold and Sevim [37 



A. Shared features 

The mechanism for selection between several interacting species is provided by the repro- 
duction rates in an individual-based population dynamics with nonoverlapping generations. 
At the end of each generation, each individual of species / reproduces asexually, giving birth 
to a fixed number F of offspring with probability Pi before dying, or it dies without off- 
spring with probability (1 — Pi). Each Pi depends on the set {nj(t)} of population sizes of 
all the species resident in the community in generation t through interspecies interactions 
and other model parameters as described below. The interactions are determined by the 



random interaction matrix M (43J, which is constructed at the beginning of a simulation 
run, and thereafter kept constant. If Mu is positive and Mji negative, then / is a predator 
and J its prey, and vice versa. If both matrix elements are positive, the relationship is a 
mutualistic one, while both negative indicate an antagonistic relationship. 

The species are defined b y a haploid, binary "genome" of length L, as in Eigen's model 
for molecular evolution [HI HH, and the potential species are identified by the index I 6 
[0, 2 L — 1]. Typically, only Af(t) <C 2 L of these potential species are actually resident in the 
community at any one time t. 

New species enter the community as each offspring organism may mutate with a small 
probability \i. Mutation consists in flipping a randomly chosen bit in the genome, and a 
mutated individual is assumed to belong to a different species than its parent, with different 
properties. Genotype and phenotype are thus in one-to-one correspondence in these models. 
This is a highly idealized picture, which is introduced to maximize the pool of different 
species available within the computational resources. The approximation is justified by a 
computational study, in which species differing by as many as L/2 bits have correlated prop- 
erties . Remarkably, it was shown that even strong correlations between the phenotypes 
of parents and offspring are relatively unimportant for the long-time dynamical properties. 

Regardless of the functional form of Pi, the time development of the mean population 
sizes, (ni(t)), is described by a set of coupled difference equations, 

(n 7 (t + l)> = (ni(t))FPi({(nj(t))m-»} 

+( f x/L)Fj2^K W (t))PK [ j ) ({(nj(t))}) , (1) 

K(I) 

where K(I) represents the L species that can be generated from / by a single mutation. 
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B. Simplified tangled- nature models 



In these models, the reproduction probability for an individual of species / is given by 
the nonlinear function 

Pl{t) = l + exp[-A 7 (i?, {nj(t)})] ' (2) 

where R is an external resource that is renewed at the same level each generation. The 
function Aj is given by 

a i e> / (+\\\ h ^ lR , \^ MijnAt) N tot{t) , . 

Here bj is a "reproduction cost" (always positive), and rjj (positive for primary producers or 
autotrophs, and zero for consumers or heterotrophs) is the ability of individuals of species / 
to utilize the external resource R, while Nq is an environmental carrying capacity [24| [a.k.a. 



Verhulst factor [461]] . The total population size is N tot (t) = Ylj n j(t)- 

Two versions of this model were studied in earlier work. In the first version there is 
no external resource or birth cost, and the off-diagonal elements of M are stochastically 
independent and uniformly distributed over [—1, +1], while the diagonal elements are zero. 
This model evolves toward mutualistic communities EL in which all species are connected 



by asymmetric, mutually positive interactions [35|, [38|, |41|, |49[ . 

In the predator-prey version of the model, a small minority of the potential species 
(typically 5%) are primary producers, while the rest are consumers. The off-diagonal part 
of the interaction matrix is antisymmetric, with the additional restriction that a producer 



cannot also prey on a consumer [34J, [35|, |37[ • In simulations we took bj and the nonzero rji as 



independent and uniformly distributed on (0, +1]. This model generates simple food webs 



with up to three trophic levels [3J, |35|, |37|. 



Both of these models provide interesting results, which include intermittent dynamics 
with power spectral densities (PSDs) of diversities and population sizes that exhibit approx- 
imate 1/ / noise, as well as power-law distributions for the lifetimes of individual species and 
the duration of quiet periods of relative evolutionary stasis. From a theoretical point of view 
they also have the great advantage that the mean-field equation for the steady-state average 
population sizes, Eq. (pQ), in the absence of mutations reduces to a set of quadratic equations 



(linear if N = oo) and thus can easily be solved exactly [35|, |37|, |38|, |49|. The models thus 
provide useful benchmarks for more realistic, but generally highly nonlinear models, like the 
ones defined below. 

The population dynamics defined by Eqs. ([2]) and ([3]) have some less realistic features. In 
particular, by summing over positive and negative terms in Aj, the models enable species 
with little food to remain near a steady state if they are also not very popular as prey, 
or have very low birth cost. Another problem is the ad-hoc nature of the normalization 
by the total population size N tot (t) in the resource and interaction terms in Aj. While 
this is the source of the models' analytic solvability, it implies an indiscriminate, universal 
competition without regard to whether or not two species directly utilize the same resources 
or share a common predator. The purpose of the present paper is to develop models with 
more realistic population dynamics and explore their complex dynamics on time scales from 
ecological (short) to evolutionary (long). 
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C. Functional-response model without adaptive foraging 



Here we develop a model with population dynamics that include competition between 
different predators that prey on the same species, as well as intraspecific competition and 
a saturation effect expected to occur for a predator with abundant prey. In doing so, we 
retain from the models discussed above the important role of the interaction matrix M, as 
well as the mutation process of the binary "genome" and the restriction to nonoverlapping 
generations. 

We first deal with the competition between predator species by defining the number 
of individuals of J that are available as prey for J, corrected for competition from other 
predator species, as 

niMu 

n " = ^pred(j) — n J > ( 4 ) 



n L M LJ 



where ^P red ( J ) runs over a n £ sucri that Mlj > 0, i.e., over all predators of J. Thus, 
^pred(j) _ .£ j . g on iy predator consuming J, then fijj = rij. 

Analogously, we define the competition-adjusted external resources available to a pro- 
ducer species I as 

R I = ^^ R . (5) 



As in the case of predators, Yli Ri = R-* an d a sole producer species has all of the external 
resources available to it: Rj = R. With these definitions, the total, competition-adjusted 
resources available for the sustenance of species / are 



prey(J) 

S I = 7 }J R J + M iJ^J 



(6) 



where ^j rey ^ 



runs over all J such that Mjj > 0, i.e., over all prey of /, and rjj = if I is 



a heterotroph. 

A central concept of the model is the functional response of species I with respect to J, 
$/j 0, 2l[. This is the rate at which an individual of species / consumes individuals of J. 



The simplest functional response corresponds to the Lotka-Volterra model |24|]: $/j = nj 



if Mu > and otherwise. However, it is reasonable to expect that the consumption rate 



should saturate in the presence of very abundant prey [2l|]. For ecosystems consisting of a 



single pair of predator and prey, or a simple chain reaching from a bottom-level producer 
through intermediate species to a top predator, the most common forms of functional re- 



sponse are due to Holling (21|. For more complicated, interconnected food webs, a number 



of functional forms have been proposed in the recent literature [§1, 0, [H, 23, 31, 32, 42], but 
there is as yet no agreement about a standard form. Here we choose the ratio-dependent 



@, B H S H [33) Holling Type II form [21], originally introduced by Getz 



1.1 



Muhu 
XSi + nj 



(7) 



where A G (0, 1] is the metabolic efficiency of converting prey biomass to predator offspring. 
The ratio dependence corresponds to intraspecific competition [la ]. Analogously, the func- 
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tional response of a producer species I toward the external resource R is 



Ad / + 

In both cases, if AS/ <C nj, then the consumption rate equals the resource {Mijhu or 
rjiRj) divided by the number of individuals of I, thus expressing intraspecific competition 
for scarce resources. In the opposite limit, XSj ^> nj, the consumption rate is proportional 
to the ratio of the specific, competition- adjusted resource to the competition-adjusted total 
available sustenance, Sj. The total consumption rate for an individual of / is therefore 

cj = * IR + P y ] $ 7J = §I = I §i/ni for A ^ <<: ni . 

1 y XSj + nj I 1/A for XSj » m 

The birth probability is assumed to be proportional to the consumption rate, 

B I = AC r J e[0,+l], (9) 

while the probability that an individual of / avoids death by predation until attempting to 
reproduce at the end of the generation is 

pred(7) 

A T = 1 - V $^ . (10) 
j 

The total reproduction probability for an individual of species / in this model is thus 

P I [t)=A I [t)B I [t) . (11) 



III. ANALYTICAL RESULTS 

The functional-response model defined in Sec. Ill CI is much less amenable to analytic 
treatment than the models we have considered previously. In particular, the simultaneous 
set of equations, 

FPj({n*j}) = 1 (12) 

where F is the fecundity, which defines the fixed-point solution {n*j} of Eq. ([fj for mul- 
tispecies communities in the mutation-free limit, cannot be solved analytically in general. 
However, some special cases can be solved explicitly. Although these analytical solutions 
are highly model-specific, they provide useful insight into some of the simplest effects of 
interspecies interactions and intra- and interspecific competition. 



A. Two competing producers with intraguild predation 

Consider two producer species characterized by their coupling constants, rji and r/ 2 . In 
the noninteracting case, M21 = Mu = 0, the species are subject to competitive exclusion, 
so that only one species, the one with the maximum value of rjj, can survive with a nonzero 
fixed-point population, n* m ^ = Ar/ max (F — 1)R. The only exception is the degenerate case 
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of 771 = 772, in which the two species are dynamically indistinguishable. The property of 
competitive exclusion in the noninteracting limit is shared with the simplified tangled-nature 
models discussed previously, and it carries over to sets of any number of noninteracting 



producers [35 



n 2 = — — u u= 1 z ^tb — mm — zi R > ( 14 ) 



For the interacting case, M 2 % £ (0, 1] (and M 12 = — M 2 i), which is a simple example of 
intraguild predation, the solution to Eq. (1121) is 

n * = ViKF - 1M(1 - M 21 ) - r|](l - M 21 ) 
1 771(1 -M 21 )[ Vl + V2 X(F-l)M 21 ]-r]i { ' 

and 

^A 2 (F-1) 2 (1-M 21 ) 2 M 21 
771(1 - M 21 )[ Vl + V2 X(F - 1)M 21 ] - 77I 

as long as both populations are nonnegative. This requires < 772 < 771 \/l — M 2 \. These 
rather complicated analytical solutions are best interpreted graphically as in Fig.[]Ja), which 
shows the case A = 1, 771 = 1, and M 2 \ = 0.5. Other parameter values give similar results. 
As r] 2 is increased from zero, the population of species 2, n 2 , first decreases weakly as 
it competes directly for resources with species 1, which is its only source of support at 
r] 2 = 0. Differentiation of the denominator in Eq. (TH|) shows that n 2 reaches its minimum 
at ?7 2 = 77^(1 — M)M/2. For larger rj 2 it increases nonlinearly due to the term quadratic in 
T) 2 in the denominator. The combined competition and predation from species 2 causes n\ 
to decrease monotonically, first linearly in r] 2 and later nonlinearly until it reaches zero at 
i] 2 = 771VI — M 2 \. For larger i] 2 , species 2 completely excludes species 1, and the stationary 
solution is n* 2 = \i] 2 (F — 1)R and n\ = 0, even though t] 2 may still be less than 771. The 
two solutions for n 2 join continuously at rj 2 = rji \/l — M 2i , and there are no other attractive 
fixed points for the mutation- free dynamics. Looking at the total population size, rii + n 2 , 
we find that it is a continuous, convex function of 77 2 , with a shallow minimum at rj 2 = 

771 |l - y/1 - X(F - 1)(1 - M)|. (The solution n\ = Xr]i(F - 1)R and n* 2 = is a fixed 

point as well, but it is repulsive under perturbations to ra 2 .) 



B. AA-species food chain 

The other case, for which the fixed-point population sizes can be found relatively easily, 
is a "food chain" in which species I + 1 feeds exclusively on the preceding one, /. The 
fixed-point equation ( fl2|) then takes the form 

F (l Ml+ ^ ) XMinU = 1 (15) 

V AM /+1 77} +n* I+1 J AM / n}_ 1 + nj ' v ; 

where for simplicity we write Mj for Mjj_i. 

With boundary conditions Uq — R and n*^ = 0, and with all Mj = M (we define 
Mi = 771), this set has a geometrically decreasing solution of the form n*j = Ra 1 with 



AM 

OL = — — 



F[AM + F(l — M) 2 } + F(l — M) — 2 j < 1 . (16) 



This solution is included in Fig. UJb) as ^ ne one corresponding to M = 00. [Parameter values 
for which a > 1 (essentially very large fecundity F combined with M and A near unity) are 
unrealistic! 
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FIG. 1: (a) Population sizes for two competing producer species with intraguild predation in 
the functional-response model without adaptive foraging, shown as functions of 772- The other 
parameters are r\\ = 1.0, M21 = 0.5, and A = 1. See the text for details, (b) Population sizes for 
food chains of M species with A = 1 and M = 0.5. See the text for details. 



If it is instead known that there are Af trophic levels, so that n*j^ +1 = 0, then the fixed- 
point equations can be solved analytically in an iterative fashion as follows. 



1. Solve the linear fixed-point equation (jlzj) expressing n*^ = const., 



for nV_i- 



2. Insert the solution for n^_ 1 in terms of rf^- into the next equation in the hierarchy 
(the one expressing n^-_ 1 = const.), 

F , ! Uv/v \, A.U A - ,// A - 



and cancel common factors to get a linear equation for n^_ 2 in terms of rfj^. 

3. Continue until obtaining Uq in terms of n*^-. 

4. Rescale the solutions to give = R. 

With large Af and /-independent Mj, this solution converges toward the decreasing geometric 
one presented above for / <C Af, as shown in Fig. [T](b). 



IV. NUMERICAL RESULTS FOR THE FUNCTIONAL-RESPONSE MODEL 

We simulated the functional-response model defined in Sec. Ill CI over 2 24 = 16 777216 
generations (plus 2 20 generations "warm-up") for the following parameters: genome length 
L — 21 (2 21 = 2 097152 potential species), external resource R = 16 000, fecundity F = 2, 
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Time (generations) Time (generations) 



FIG. 2: Time series of diversities (a) and population sizes (b) for one specific simulation run 
of the functional-response model without adaptive foraging. The strongly fluctuating curves in 
the background are sampled every 8192 generations, while the smooth curves with data points in 
contrasting colors that are overlaid in the foreground are running averages over 524 288 generations. 
Black with light gray (yellow online) overlay: all species. Light gray with dark gray overlay 
(green and magenta online): producers. Dark gray with light gray overlay (red and cyan online): 
consumers. 

mutation rate fi = 10~ 3 , proportion of producers c pro d = 0.05, interaction matrix M with 
connectance C — 0.1 and nonzero elements with a symmetric, triangular distribution over 
[—1,-1-1], and A = 1.0. (The high value of A is of course biologically unrealistic, and it 
was chosen to obtain a larger population of heterotrophs for a computationally manageable 
autotroph population.) We ran five independent runs, each starting from 100 individuals of 
a single, randomly chosen producer species. 



A. Time series 



Time series of diversities (effective numbers of species) and population sizes for one 
realization are shown in Fig. [2j To filter out noise from low-population, unsuccessful 
mutations, we define the diversity as the exponential Shannon- Wiener index j20|. This 
is the exponential function of the information-theoretical entropy of the population dis- 
tributions, D(t) = exp[S({ni(i)})], where S ({nj(t)}) = - J2{i\ PI (t)>o} Pitt) m M*) witn 
Pi(t) = nj(t)/N tot (i) for the case of all species, and analogously for the producers and 
consumers separately. 

The time series for both diversities and population sizes display intermittent behavior 
with quiet periods of varying lengths, separated by brief periods of high evolutionary activity. 
The intermittency is highlighted by the time series for the accumulation of new and extinct 
species, shown in Fig. [3j In this respect, the results are similar to those seen for the simplified 
tangled-nature models in earlier work [34]) 35], H, 41]. However, diverse communities in this 



model are less stable than those produced by the simplified tangled-nature models. It 
is possible that this instability is related to the tendency of "triangles" consisting of two 
species competing for a common resource while one of them also feeds on the other one 
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FIG. 3: (a) Time series of the accumulation of new and extinct species in the same simulation 
run depicted in Fig. [2j The solid, black curve shows the total number of different species that 
have at least once attained a population size nj > 1000 by time t. The dashed curves count the 
total number of species that have gone extinct after attaining a maximum population greater than 
1000. The black dashed curve refers to all species, the light gray one (green online) to producers, 
and the dark gray one (red online) to consumers. The ratio of approximately 1.89 between the 
dashed and full black curves indicate that major species recur on average about twice during the 
evolution. This is an artifact of the finite genome length, (b) The detailed, intermittent structure 
of new species (the solid, black curve in (a)) over 400 000 generations. The interval is indicated by 
the short, horizontal bar in (a). 

(intraguild predation) to collapse, that we discussed in Sec. IIII Al The instability expresses 
itself in a tendency for this model to flip randomly between an active phase with a diversity 
near ten, and a "garden of Eden" phase of one or a few producer species with a very low 
population size of numerous unstable consumer species, such as the one seen around 10 
million generations in Figs. [2]and[2](a). 



B. Power-spectral densities 

To obtain information about the intensity of fluctuations in the evolving community, we 
calculate power-spectral densities, or PSDs.|50j] These are presented in Fig. H]for the diver- 
sities and the population sizes (Fig. BJa)) and the intensity of extinction events (Fig. IH(b)). 
The former two are shown for the total population, as well as separately for the produc- 
ers and consumers. All three are similar. Extinction events are recorded as the number of 
species that have attained a population size greater than one, which go extinct in generation 
t (marked as "Species" in the figure), while extinction sizes are calculated by adding the 
maximum population sizes attained by all species that go extinct in generation t (marked as 
"Population" in the figure). The PSDs for all the quantities shown exhibit approximate 1/f 
behavior. For the diversities and population sizes, this power law extends over more than 
five decades in time. The extinction measures, on the other hand, have a large background 
of white noise for frequencies above 10~ 3 generations -1 , probably due to the high rate of 
extinction of unsuccessful mutants. For lower frequencies, however, the behavior is consis- 
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Frequency (1/generation) Frequency (1 /generation) 



FIG. 4: (a) PSDs for the diversities and population sizes, each recorded separately for all species 
and for producers and consumers. The time series were sampled every 16 generations, (b) PSDs 
for the extinction activity, sampled every generation. In both parts of the figure, the results are 
averaged over five independent simulation runs. See discussion in the text. 



tent with 1/f noise within the limited accuracy of our results. We note that the apparent 
1/ / behavior in the PSD of extinction events in the model of Nunes Amaral and Meyer [28[ 
extends over only one decade in frequency. 



C. Species lifetimes and durations of quiet periods 



The evolutionary dynamics can also be characterized by histograms of characteristic time 
intervals, such as the time from emergence till extinction of a species (species lifetimes) or 
the time intervals during which some indicator of evolutionary activity remains continu- 
ously below a chosen cutoff (duration of evolutionarily quiet periods). Histograms of species 
lifetimes are shown in Fig. 0(a). As our indicator of evolutionary activity we use the magni- 
tude of the logarithmic derivative of the diversity, |di?/dt|, and histograms for the resulting 
durations of quiet periods, calculated with different cutoffs, are shown in Fig. 0(b). Both 
quantities display approximate power-law behavior with an exponent near —2, consistent 
with the 1/f behavior observed in the PSDs 



29, 38 



It is interesting to note that the distri- 
butions for these two quantities for this model have approximately the _same exponent. This 
is consistent with the previously studied, mutualistic model & 



38 , but not with the 



predator-prey model [34], [35|, [37j. We believe the linking of the power laws for the species 
lifetimes and the duration of quiet periods indicate that the communities formed by the 
present model are relatively fragile, so that all member species tend to go extinct together 
in an avalanche-like "mass extinction." In contrast, the previously studied predator-prey 
model produces simple food webs that are much more resilient against the loss of a few 
species, and as a result the distribution of quiet-period durations decays with an exponent 



near 



3UI22I. 
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FIG. 5: (a) Histograms of species lifetimes, shown for all species, as well as separately for producers 
and consumers, (b) Histograms of the durations of evolutionarily quiet periods, defined as the 
times that the logarithmic derivative of the diversity, IdS'/dtj (averaged over 16 generations), falls 
continuously below some cutoff. The inset is a histogram of dS/dt, showing a Gaussian center with 
approximately exponential wings. The parabola in the foreground is a Gaussian fit to this central 
peak. The cutoff values for the main figure, between 0.008 and 0.024, were chosen on the basis of 
this distribution. The data in both parts of the figure are averaged over five independent simulation 
runs, and the error bars represent standard errors based on the differences between runs. 
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FIG. 6: Population sizes for two competing producer species with intraguild predation in the 
model with adaptive foraging, shown as functions of rj2- (a) 771 = 1.0. (b) 771 = 0.5. In both cases, 
A = 1 and M21 = 0.5. See further discussion in the text. 



V. MODEL WITH ADAPTIVE FORAGING 

The model studied above is one in which species forage indiscriminately over all available 
resources, with the output only limited by competition. Also, there is an implication that an 
individual's total foraging effort increases proportionally with the number of species to which 
it is connected by a positive M u . A more realistic picture would be that an individual's 
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total foraging effort is constant and can either be divided equally, or concentrated on richer 
resources. The latter constitutes adaptive f orag ing. While one can go to considerable length 
devising optimal foraging strategies 0, 0, fell 32J, we here only use a simple scheme, in 



which individuals of / show a preference for prey species J, based on the interactions and 
population sizes (uncorrected for interspecific competition). The proportion of its foraging 
effort that / allots to J is approximated as 



M u nj 
ViR + Ex ey(/) M IK n K 
and analogously for the effort assigned to the external resource 

ViR 

ViR + E P K 7il> M IK n K 



Pu = ~ , _ prcy(J) ; : , (17) 



Pir = - _ p / cy(J) ; ; ■ (18) 



The total foraging effort is thus normalized: Pm + Yl^ 7 ^ Pu = 1- These preference factors 
are used to modify the reproduction probabilities by replacing all occurrences of Mjj by 
Mijpu and of 77/ by T]ip IR in Eqs. (SHE]). 

The adaptive foraging obviously has no effect on a simple food chain since no species in 
this case has more than one choice of prey. The analytical results thus remain as discussed 
in Sec. IfflBl 

For the case of two competing producers with intraguild predation we did not obtain 
analytical results for the fixed-point population sizes, except for the special cases of r\ 2 = 0, 
which corresponds to a two-species food chain, and of M 2 \ = 0, which reduces to simple 
competitive exclusion of the species with the lower t]j. However, numerical results for M 2 ± = 
0.5, obtained by iteration of Eq. (CQ) with /1 = 0, are given in Fig. [6j The results are 
similar for other values of the model parameters. The parameters in Fig. El^a) are the 
same as in Fig. Q](a), and we see that the regime of two-species coexistence is significantly 
extended by the adaptive foraging and here covers the full range of r/ 2 . We also see that 
while n 2 is reduced, compared to the case without adaptive foraging, both n\ and the total 
population, n\ + n 2 , are significantly increased, indicating a more efficient overall resource 
utilization by the community. In Fig. Mjo) we reduce r\\ to 0.5 to explore the possibility that 
772 > 771. We find that the coexistence solution extends up to r\ 2 ~ 0.553, where n\ = n 2 
and the solution changes discontinuously to the familiar n\ = \(F — l)r] 2 R and n\ = 0. 
In fact, for rj 2 between r\\ = 0.5 and 0.553, both solutions are locally stable under small 
perturbations. This is indicated by the dotted lines in the figure. Outside this range, the 
solutions shown are globally stable attractors. As in the case without adaptive foraging, the 
solution n* = \r)x(F — 1)R and n 2 = is repulsive under perturbations to n 2 . 

The results of implementing the adaptive foraging strategy in long-time simulations of 
evolving multispecies communities are quite striking. The system now has a metastable low- 
diversity phase similar to the active phase of the non-adaptive model, from which it switches 
at a random time to a stable high-diversity phase with much smaller fluctuations. As seen in 
Fig. El the switchover is quite abrupt, and Fig. [S] shows that it is accompanied by a sudden 
reduction in the rate of emergence of new species. The existence of a stable high- diversity 
phase with increased producer and total populations and reduced consumer population in 
the case of adaptive foraging is consistent with the extension of the stability of two-species 
coexistence discussed above. The increased total population is consistent with the improved 
resource utilization, and the sudden nature of the transition from low to high diversity is 
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FIG. 7: Time series of diversities (a) and population sizes (b) for the model with adaptive foraging. 
The interpretation of the colors and lines are the same as in Fig. [2j 
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FIG. 8: (a) Time series of the number of new species that have reached a population greater 
than 1000 (lower curve) and greater than 100 (upper curve) in the same simulation run depicted 
in Fig. [71 (b) The detailed, intermittent structure of the upper curve in (a) on a very fine scale of 
6000 generations. See discussion in the text. 



consistent with the discontinuity and bistability, all observed in the adaptive two-species 
case. As adaptive foraging implies an effective reduction of omnivory, our results present a 
scenario in which reduced omnivory leads to increased community stability [251 ] . 

As seen in Fig. [9j the PSDs for both the diversities and population sizes in both phases 
show approximate 1/f noise for frequencies above 10~ 5 generations -1 . For lower frequencies, 
the metastable phase shows no discernible frequency dependence, while for the stable phase, 
the frequency dependence continues at least another decade. It thus appears that long-time 
correlations are not seen beyond about 10 5 generations in the metastable phase. 

In fact, the system can also escape from the low-diversity phase to total extinction, which 
is an absorbing state, and in some of our simulation runs we avoided this by limiting \Mjj\ 
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FIG. 9: PSDs for the diversities and population sizes in the metastable phase (averaged over 
three independent runs) (a), and the stable phase (averaged over five independent runs) (b) for 
the model with adaptive foraging. Both show approximate 1/f noise for frequencies above about 
10 -5 . However, in the metastable phase the PSDs appear to approach constant levels for the lowest 
frequencies, indicating the absence of long-time correlations. 



to less than 0.9. This restriction does not seem to have any effect on the dynamics in the 
high-diversity phase. 



VI. CONCLUSIONS 



In this paper we have extended our study of the long-time dynamics of a class of 
individual-based models with stochastic population dynamics, nonoverlapping generations, 
and random mutations during reproduction. Previous studies concentrated on simplified 
versions of the tangled-nature model 0, ?l 17], sharing with that model a reproduction 



probability denned by Eqs. Q2p and © [34], |35|, |37|, |38|, |41|, |49j. While these models in 
the absence of mutations allow for exact, analytical solutions for the fixed-point popula- 
tions of any given community of species, this convenient mathematical property is due to 
two somewhat unrealistic features, (i) The lumping together of resources and predation 
as respectively positive and negative contributions to the quantity Aj in Eq. (ii) The 
normalization of the population sizes nj in Aj by the total population size A^ to t, which can 
be seen as an indiscriminate, universal competition effect. 

Here, we therefore introduce population dynamics based on a functional response for 
predators versus their prey and for autotrophs versus the external resource. Interspecific 
competition is introduced through the competition-adjusted resources defined in Eq. 
and intraspecific competition is accounted for by the ratio-dependent functional response 
due to Getz [l6|, Eq. ([7]). The probability for a live individual to give birth is given by 
Eq. (Q, and the probability that an individual avoids being eaten before it can reproduce 
is given by Eq. ffTUj) . Their product, Eq. ffTTj) . is the total reproduction probability Pi. 

We considered two versions of this model, one without adaptive foraging, and one with 
it. While a complete analytic solution of the fixed-point communities is not feasible for 
either model, we obtained solutions for a simple food chain of predators supplied by a single 
producer species. For the model without adaptive foraging we also obtained an analytic 
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solution for the coexistence of two producer species, one of which also acts as a predator 
toward the other (intraguild predation). A corresponding analytic solution was not found 
for the model with adaptive foraging. However, numerical solutions of Eq. ([1]) with [i = 
showed a significantly expanded parameter range for coexistence, including a regime where 
both coexistence and competitive exclusion are locally attractive fixed-point solutions. The 
coexistence solution also exhibited a decrease in the predator population, which was more 
than compensated by an increase in the prey population, leading to a significant increase 
in the total population size; in other words to a more efficient resource exploitation by the 
overall community. 

Long-time kinetic Monte Carlo simulations of the two models produced time series of 
diversities and population sizes that exhibit approximate 1/f noise in a wide frequency 
range. However, there are significant differences between the models. 

Without adaptive foraging, the community flips randomly back and forth between an 
evolutionarily active phase with diversities around ten, and another phase with a very low 
evolutionary turnover, a small number of coexisting producer species, and a very small 
population of unsuccessful consumers. 

With adaptive foraging, the model displays a metastable phase which resembles the active 
phase in the previous case. After a random amount of time, this phase suddenly gives way, 
either to total extinction, or to a new, stable phase with an order of magnitude larger 
diversities and somewhat higher total population size. The relative fluctuations in this 
phase are much smaller than in the other phases in either model. We find it reasonable to 
believe that the increased diversity and population size in the stable phase of the model 
with adaptive foraging are related to the increased tendency toward species coexistence in 
this model, observed in Fig. EJ This stabilization of the communities is consistent with 
observations for the web-world model with adaptive foraging [si, 0, 31, 32]. The dynamics 



and community structures of the stable phase in the model with adaptive foraging are studied 
in detail in a forthcoming paper 36| . 



In summary, the models studied here show that some aspects of the long-time dynam- 
ics, such as the 1/f noise in diversities, population sizes, and extinction sizes, are quite 
robust. However, the community structures and their stability or lack thereof show signif- 
icant differences, both from previously studied tangled-nature type models, and depending 
on whether or not adaptive foraging is implemented. A comprehensive understanding of 
the universal and nonuniversal fluctuation properties of large-scale coevolution and their 
relation to avalanches of species extinctions is likely to demand a combination of implemen- 
tations of more realistic population dynamics and mutation mechanisms, and further study 
of minimalistic models. 
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